Weak ferromagnetism and spiral spin structures in 
honeycomb Hubbard planes 



MAN Araujo 1 2 and N M R Peres 2 - 3 

1 Departamento de Ffsica, Universidade de Evora, P-7000-671, Evora, Portugal 
2 Center of Physics, Universidade do Minho, P-4710-057, Braga, Portugal 
3 Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, MA 
02215,USA 

E-mail: mana@uevora.pt 

Abstract. Within the Hartree Fock- RPA analysis, we derive the spin wave spectrum 
for the weak ferromagnetic phase of the Hubbard model on the honeycomb lattice. 
Assuming a uniform magnetization, the polar (optical) and acoustic branches of the 
spin wave excitations are determined. The bipartite lattice geometry produces a q- 
dependent phase difference between the spin wave amplitudes on the two sub-lattices. 
We also find an instability of the uniform weakly magnetized configuration to a weak 
antiferromagnetic spiraling spin structure, in the lattice plane, with wave vector Q 
along the T — K direction, for electron densities n > 0.6. We discuss the effect 
of diagonal disorder on both the creation of electron bound states, enhancement of 
the density of states, and the possible relevance of these effects to disorder induced 
ferromagnetism, as observed in proton irradiated graphite. 
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1. Introduction 

Recent interest in strongly correlated systems in non-square lattices, such as the 
triangular, honeycomb and kagome lattices, is justified by the possible realization of 
exotic metallic [HE], magnetic Ej and superconducting states [Hj both in inorganic 
and organic materials. From the organic side, graphite and related carbon allotropes 
are physical systems where growing evidence for exotic types of ground states is being 
accumulated during the last few years. In graphite, for example, experimental research 
put forward evidence for unusual metallic and magnetic properties El El El- In 
particular, ferromagnetism has been observed at high temperature in graphite |Hj 
which may not be due to magnetic impurities. Also, the observation [2] of magnetic 
order induced by proton irradiation challenges the theoretical description. Graphite is 
not alone on the ferromagnetic-order-by-disorder scenario, with the inorganic CaRuO"3 
material also exhibiting disorder-induced ferromagnetism [S] . Recent experimental work 
[HI EH E] has produced atomic thin graphite planes where the exciting physics of 2D 
Dirac fermions may be directly observable. Motivated by these experimental studies 
and because the microscopic origin of ferromagnetism in these compounds is far from 
being understood, we decided to study the magnetic properties of a doped Hubbard 
model on an honeycomb lattice - a single graphite plane. To the best of our knowledge, 
ferromagnetic spin waves in the honeycomb lattice (as an itinerant electron system) 
have not been studied in the past. The fact that the honeycomb lattice is a Bravais 
lattice with a basis immediately presents us the possibility of observing both polar 
and acoustic spin waves ^21 EI]- We focus our research on the stability of the weak 
homogeneous ferromagnetic phase found in [14] . The paper is organized as follows: the 
model is introduced in section and the energy spectrum for the weak homogeneous 
ferromagnetic system is derived; Section El is devoted to the description of spin waves 
in the weak homogeneous ferromagnetic phase; the possibility of spiral spin states is 
investigated in Section 01 where a spiral arrangement is found with lower energy than 
that of a uniform magnetization; Section El gives a discussion on the possibility of 
formation of electronic bound states, due to impurities, and of the possible relevance of 
disorder to experiments on proton irradiated graphite. 

2. Model Hamiltonian 

The Hubbard model is defined as 



where c\ a (q )(7 ) represents a creation (destruction) electron operator with spin a at site i, 
tij is the hopping integral between two sites i and j, U is the on-site Coulomb repulsion, 
and fi is the chemical potential. In the honeycomb lattice we identify two sub-lattices, 
A and B (see figure Q), where the primitive vectors of the underlying triangular lattice 
are denoted by a\ and a 2 . For later use we also define the vector = (cti — a 2 ). The 
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reciprocal lattice vectors are 61 and b 2 and define a hexagonal shaped first Brillouin 
Zone. Also shown in figure ^ are the vectors connecting any A atom to its nearest 
neighbours, denoted as S\, 62 and #3. The electrons leaving on each sub-lattice will be 




Figure 1. Primitive vectors (ai and 0,2) for the honeycomb lattice. The vectors 6±, 
82 and <$3 connect the A site to its three neighbouring B sites. The hexagonal first 
Brillouin zone corresponds to the reciprocal lattice vectors, b\ and hi- 



denoted by field operators a and 6, respectively. The Fourier transformation between 
real and momentum spaces is given by: 



= E ^4,. > b l = 4r E e ih - (2) 



a f = 

J.O" /AT / J ~ —K,(T ' ~1,<T /AT 



and we take Na = Nb = iV as the number of unit cells. In the calculations below, we 
shall consider first and second neighbour hopping integrals, t and t', respectively. The 
Hubbard model then takes the form: 



H = EP( fc ) " f*)(4,„ak,* + bljk,*) 



fc,cr 



+ E[<K fc )4,A,. + 0*(*OC a ^] + H v , (3) 



fc,cr 



with 



D(k) = - 2t'J2 C0S ( a i ■ k ) 



i=l 

3 



m = -*$>^. (4) 



i=l 

In the ferromagnetic ground state, the average occupancy of lattices sites is given by 
. + n m /lt , . n m , N 

«aOi,<r) = 2 + ^ y ' (%A«r) = g + a y ' ( 5 ) 

with the spin index a = ±1. This may also be generalized to describe antiferromagnetic 
ordering if we replace m with — m in one of the equations © 14j . An Hartree-Fock 
treatment of the Hubbard term, Hu, taking into account equation ©, yields a set of 
quasi-particle bands given by 

E%(k) = D(k) + ^(n - am) + a\<f> h \ , (6) 
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where a = ± is a band index. In the ferromagnetic phase the single particle Green's 
functions can be written, in momentum space, as: 

g?(iu n ,k)= D . , (7) 

j= ± %u) n - Eij{k) 

A p iS{k) In 

g b a a (iuj n ,k) = v — p- (9) 

j= ± iu n - Ei(k) 

gl b (iuj n ,k) =g a a a (u n ,k), (io) 

where we have defined e l5 ^ = 0(fe)/|0(fc)|. The Hartree-Fock magnetization is given 
bym=(l/2N)E^(Tf[E^k)]. 

Since we are mainly concerned with the ferromagnetic phase, the calculations 
in sections and H] are performed for an electronic density smaller than half filling. 
Therefore, electrons at the Fermi level will not be treated as massless Dirac fermions. 
Such treatment is usually appropriate for electrons in graphite planes at half filling (or 
close to half filling) [TH] . 

3. Magnetic collective excitations 

We obtain the magnetic collective excitations from the poles of the transverse spin 
susceptibility calculated in the RPA approximation. Because there are two sub-lattices, 
the susceptibility is actually a second order tensor given by the expression 

&.(q,iu> n ) = fJ T ' dTe**> r (T T S?(q,T)Sr(-q,0)) (11) 

where i,j = a,b are sub-lattice labels and Sf(q), S~(q) are the spin-raising and lowering 
operators for each sub-lattice. The RPA expansion gives a Dyson equation for the 
transverse spin susceptibility, which can by written, in matrix form, as 



x{q, tun 



u 



l 



i--x\q^n)\ x^q^n) (12) 

where 1 denotes the 2x2 identity matrix. The poles of the susceptibility tensor, 
corresponding to the magnetic excitations, are then obtained from the condition: 

det[l-^x°(q,oo + i0 + )} =0. (13) 

Below the particle-hole continuum of excitations, the spectral (delta-function 
contributions) part in x+- ''{Qi^ + vanish and there is the additional relation 
X+- — (x+- a6 )*- The zero order susceptibility tensor, x+_(qr, u>), can be written as 

4 k-,a^ 2 =± iUJ n + E ] (k) -E^ik-q) 
where the matrix A(k, q) contains the coherence factors: 
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(3(k,q) = exp[— i5(k) + i5(k — q)], with f" 2 (k) representing the Fermi function with 
argument E^ 2 (k), and equivalent representations hold for the other cases. 

In addition to the collective ferromagnetic spin waves there are also single particle 
flip-spin excitations which define the so called Stoner continuum. In our case we may 
have up to four regions in the energy-momentum plane associated with the latter type 
of excitations. Their spectra are given by 

A ai ' Q2 (q) = E^(k -q)- E* 2 {k) . (16) 

Depending on the position of the Fermi level, one or two of the regions defined by 
equation (|16|) may not occur because the bands may be empty. Because we have 
the lower (E~) bands partially filled, there always are two Stoner regions given by 
A~~(q) and A + ~(q). We consider only this case below (for the other cases the results 
are qualitatively the same). In figures El and El the left panel shows the the Stoner 
continuum defined by A~~(q) as the area enclosed by the dashed-dotted line starting 
at A~'~(0) = Um. The model parameters have been chosen so that the system is 
definitely not antiferromagnetic. In figures El and El we plot the solutions to equation 
(|13|). It is clear that the system has two different types of collective magnetic excitations, 
namely the usual acoustic mode u) ac {k) and the polar or optical mode u opt (k), associated 
with the existence of two atoms per unit cell ^21 fTTTj . It is quite interesting that the 
behaviour of both branches of excitations does not follow the same trend along different 
directions of the Brillouin zone: along the T — X direction (see figure EJ), for example, 
and for \q\ > 0.25rX only the optical branch remains. On the other hand, in the 
r — K direction (see figure EJ) it is the optical branch that vanishes at \q\ ~ 0.2TK 
while the acoustic mode survives. The vanishing frequency of the acoustic (or optical) 
modes at finite momentum is associated with an instability of the homogeneous weak 
ferromagnetic phase toward a state exhibiting possibly weak ferromagnetic order in the 
z direction and spiral order in the xy plane, which will be analyzed in section 0] below. 

The eigenvector of the matrix x^(Qi ^(q)) that is associated with the eigenvalue 
N/U gives the spin wave amplitudes over the A and B sub-lattices. We note that 
equation (JT5j) is equivalent to the eigenproblem 



" (si) ' 


1 


' (si) ' 


. (s^) _ 


" u 


. (si) . 



where (Si), (S^) denote the spin wave amplitudes over the two sub-lattices [TH]. In our 
case, ^(°)' aa = and x^ ab = (x^ ba )*- The equations for the eigenvalue A and the 

corresponding eigenvector are: 

A = x (0) ' aa ± lx (0)a6 | (18) 

(O)afe 

(^> =± §^|<^>. ( 19 ) 

where A = N/U is the relevant eigenvalue, as can be seen from equation (fTTj) . Equation 
()19|) shows that the spin wave amplitudes are related by a phase factor. Therefore, the 
phase of the complex matrix element x^ ab determines the angle between the transverse 
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Figure 2. Left panel: Spin waves (line with circles for u> op t and continuous line for 
Lo ac ) and Stoner continuum border (thin dashed-dotted line) along the T — X direction 
of the Brillouin zone. Right: band energies along the T—X direction (the band energies 
are measured relatively to the chemical potential); The parameters are: U = 4, t — 1, 
t' = —0.2, n = 0.75, and the magnetization and chemical potential are m = 0.25 and 
/i = 0.36. The energy through all this paper is in units of t. (Subscripts u and d in the 
right 




Figure 3. Left panel: Spin waves (line with circles for uj opt and continuous line for 
Lo ac ) and Stoner continuum border (thin dashed-dotted line) along the T — K direction 
of the Brillouin zone. Right panel: quasi-particle band structure along the V — K 
direction (band energies measured relative to chemical potential) ; The parameters are 
the same as in figure [21 (Subscripts u and d in the right panel denote spin projections 
t and I, respectively.) 



spin components (Sj[) and (<S#). The optical mode in the T — X direction (shown in 
figure EJ) starts off with (Sj[) = —(Sr) for small q, as expected of an optical mode, 
but the angle between (S^) and (S^) monotonically decreases from tc, upon increasing 
wave vector, and equals n/2 when to opt attains its minimum. At that point, x^ ab is 
pure imaginary. On the other hand, in the acoustic mode the angle increases from zero 
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Figure 4. Relative position of the precessing spins on neighbouring A and B sites for 
a spin wave in the T — X direction: (a) acoustic mode with small q; (b) optical mode 
with small q; (c) minimum (maximum) frequency of the optical (acoustic) mode. 

to 7r/2, when u ac is maximum (the angle between the spins is illustrated in figure H}. 
As the wave vector further increases, the acoustic mode frequency rapidly decreases and 
vanishes at q « 0.25rX. 

Considering now the spin waves in the V — K direction, the angle between the 
precessing spins is always zero or ir in the acoustic and optical modes, respectively. The 
optical mode frequency vanishes shortly after the interception with the acoustic branch 
and only the latter survives for increasing q. 

The phase difference between spin wave amplitudes just described for the T — X 
direction is a manifestation of the complex coherence factors appearing in the single 
particle Green functions (JHJ) and Q, resulting from the honeycomb lattice geometry. 

4. Spiral spin states 

The disappearance of the acoustic mode, at wave vector q « 0.25TX, and of the optical 
mode, at wave vector q « 0.20r.K~, suggests an instability to a spiral spin state [T7j . 
In such a state, the spiral spin configuration is characterized by non-zero transverse 
magnetization at site j, (Sj) = (S + )e l( ^' Rj , in addition to a uniform alignment in the z 
direction, (S z ). The amplitudes of the spiral, (Sjum), at sub-lattices A and B are given 
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by: 



(St) 



TT J2( a k,-] a k+Q,i) 
JV C k 



(20) 



In general, there will be a nonzero angle # between the transverse sub-lattice 
magnetizations (S^) and (S B ), so that (5*^) = e l9 {S\). 

Following the mean field equations are obtained from the minimization of 
the ground state energy with respect to the order parameters (S\^) and (S z ). Each 
Bloch fe-state, 7*., representing an elementary excitation, is a linear superposition of 
the fields a k ,i, afc+Q,j.; ^fc,T an d ^fc+Q,|- Conversely, we can rewrite each of the fields as a 
combination of Bloch states, and recast the expectation value of the kinetic term in Q 
as well as the order parameters ()20|) in terms of 7*. operators. Using Wick's theorem, 
the expectation value of the Hubbard term in the Hamiltonian (jHJ) can be expressed as: 

(21) 



(Hu) = U(j - (S Z A ) 2 - (S z B y - (S2){3t) ~ (Sb){S£)). 



By minimizing (H), as given in equations and (|21|). we find that the Bloch fc-state, 
7^, diagonalizes an effective 4x4 Hamiltonian matrix, H e ff(k), which can be expressed 
in the basis (a fcjT , a k+Qil , b kth b k+Qii ), as 

V A (k) C*(k) 



H ^ {k) L C{k) V B (k) 
where the 2x2 matrices D a (with a — A, B) and C are given by: 

~ D(k)-U(S»)-p -U(S+) 



V a (k) 



-u(s-) 



D{k + Q) + U(S z ) - n 



(22) 



(23) 



and 



C(k) 



0(fc) 




(24) 





(fc + Q) 

The mean field equations only determine the phase difference between (S^) and (S B ), 
so we choose (Sj[) to be real. At each point in the weak ferromagnetic region of the 
phase diagram we have to choose the spiral wave vector Q that minimizes the ground 
state energy. This vector should lie along one of the high symmetry directions in the 
Brillouin Zone. 

We first look for the most favorable wave vectors lying along the T—X direction. For 
the same parameters as in figure El we find a spiral state with Q = \TX , (S z ) = 0.12 
and (S^) = 0.036. We note that this spiraling state has a smaller ^-component of 
the magnetization than that in the uniform phase. The angle between transverse 
magnetizations 9 = 0.77n ~ 37r/4. We find, however, that the energy difference between 
this spiral state and that with uniform magnetization is indeed very small, not exceeding 
10 -4 £ per lattice site. This has been checked for several lattice sizes. We also note that 
the obtained spin transverse component is not negligible compared to (S z ). The most 



Weak ferromagnetism and spiral spin structures in honeycomb Hubbard planes 



9 



favorable spiral wave vector depends slightly on interaction and density: at U = 3.5 and 
n = 0.8, for instance, Q = §TX is the most favorable, with (S z ) = 0.089, (S%) = 0.046 
and 9 = 0.73tt. 

We now consider states with Q oc TK. Overall, the energies of these spiral states 
are found to be lower than those of the states with Q oc TX considered above. For the 
same parameters as in figure El we find the optimum wave vector to be Q = f TK, where 
(S z ) vanishes and (Sj[) = 0.24 = —(S^). The spin configuration is, therefore, planar 
and the two sub-lattices have opposite magnetizations. The energy (per lattice site) 
is 0.02 lower than that with uniform magnetization. An approximate representation of 
this state is shown in figure El The length of the most energetically favorable Q depends 




Figure 5. Representation of the spin transverse components (5| B ) for the spiral 
states with: Q = (1/4)TX (left); Q = (3/i)TK (right). The latter configuration is 
purely planar and has the lowest energy for the same parameters as in figure 

significantly on U and n. If, for instance, U = 3.5 and n = 0.8, then Q = |TK is the 
most favorable, with (S^) = 0.18 = — (S^). The energy per lattice site of this state is 
0.008 lower than that of the uniform state. In such a planar spin configuration, a small 
ferromagnetic alignment along z could still arise in the presence of a weak anisotropy 
or external magnetic field. 

An anisotropic perturbation producing an easy axis (or a small magnetic field) need 
not exceed an energy of the order 10~ 2 t per lattice site in order to induce the uniform 
state, specially at the lower U values considered. It is still possible that the minimum 
energy spin structure of the system could be a superposition of spirals with different Q 
vectors, as has been found for the n = 1 Hubbard triangular lattice jTHj. The search for 
such structures, as well as their sensitivity to disorder, is beyond the scope of this work, 
however. We also find that the spiral states are absent at smaller densities (n < 0.6). 
A schematic representations of our findings is shown in figure EH 

5. Disorder as possible mechanism to ferromagnetism 

Ferromagnetic order induced by disorder in proton irradiated graphite was observed by 
Y. Kopelevich et al. [2j. The disorder induced by proton irradiation can be modeled 
by diagonal disorder, where the local energy of some sites is modified. The problem of 
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Figure 6. Schematic phase diagram of a honeycomb layer, with t' ^ 0, showing 
the paramagnetic (P), antiferromagnetic (AF), weak uniform ferromagnetic (WF) and 
spiral phases. The Nagaoka phase is a fully polarized ferromagnet. The exact position 
of the transition lines is not given, except at some special values marked in the axes. 
The results for t' = are not qualitatively different. 



treating disorder and Coulomb interaction together is a hard one in condensed matter 
physics [19]. Therefore we start by studying the effect of a finite density of uncorrelated 
impurities on the electron gas in the honeycomb lattice at half filling, where the low- 
energy electronic excitations can be described by massless Dirac fermions ^Hj- The 
effect of disorder on the properties of Dirac fermions leads to some unexpected results, 
and was discussed in the context of disordered superconductors by some authors [23 EE] 
In our study, we compute the T— matrix using the massless Dirac fermion description. 

The Matsubara Green's functions are determined via the equation-of-motion 
method. After the usual manipulations |22| we obtain the 2x2 Green's function matrix 
as 

G(p, iuJn) = G°(p, iu n ) + G°(p, iuj n )T{iuj n )G°{ P) iuj n ) , (25) 
with the Matsubara T— matrix given by 



where G°(p,iuj n ) is the Green's function matrix with V = 0, and 
G°{iu n ) = ^G°(p,iu; n ). 



(26) 



(27) 



For a small density, n imp , of scatterers (but finite in the thermodynamic limit), and 
after the position of the scatterers has been averaged over ensemble configurations, the 
Green's function matrix can be written as 

G(p,iu n ) = [[GP{p,iu n )\- x - E(zuv)]- 1 , (28) 
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where |221 

S(iw n ) = Vn imp [l - ^G\iu n )}- 1 . (29) 

Once the Green's function (|28|) is known, one can proceed to include the effect of 
correlations into the problem. The electronic bound states are given by the poles of 
(J2HJ). There always is a bound state due to the impurity, below the energy band, 
independently of the value of V . 

The existence of bound states allows for a possible mechanism to the disorder 
induced ferromagnetic behaviour in proton irradiated graphite. Graphite is usually 
modeled as a half-filled honeycomb plane, where electrons near the Fermi level have 
linear dispersion ^S] • The sample irradiation produces the displacement of the carbon 
atoms from their original position. In this case, even if hydrogen atoms become bonded 
to some of the carbons, from the lattice point of view a dilution of lattice points is being 
induced. In this case, we should take the limit V — ► oo, even if rzj mp is small. This effect 
leads to drastic change in the density of states p(w), where a strong enhancement of 
p{uS) in the vicinity of u = is obtained. Such an enhancement can be responsible for a 
large reduction of the critical U needed for ferromagnetism, as follows from the Stoner 
criterion. The effect of disorder on the density of states of Dirac fermions is shown in 
figure [7| for several values of n imp and V. If V is negative there are bound states for 
the electrons below the bottom of the band. As V increases an enhancement of p(u) in 
the vicinity of u = starts to develop. We have, therefore, two possible routes toward 

0.4 
0.3 
§ 0.2 

Q. 

0.1 



"-3 -2 -1 1 2 3 -0.2 0.2 

w 0) 

Figure 7. Density of states p(u>) for several impurity densities rii mp and V values, in 
a model of Dirac fermions. The left panel show p over the bandwidth. The right panel 
shows a zoom of the density of states presented in the left panel. 

the appearance of magnetic order in graphite, namely bound states and enhancement of 
the density of states. This type of enhancement of p(u) is characteristic of acoustic 
excitations, either fermionic (Dirac fermions) or bosonic fmagnons) |23| . This very 
qualitative view has to be corroborated by more detailed calculations taking into account 
both disorder and Coulomb interaction. In particular, it is important to compute how 
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the critical lines from the paramagnetic to the magnetic phases change with the amount 
of disorder. These issues will be the subject of a future publication. 

6. Conclusion 

In conclusion, we have studied the spin collective excitations of the homogeneous weak 
ferromagnetic state in the honeycomb lattice and found an instability to spiral spin 
structures at electron densities above n « 0.6 Although our calculation is performed 
for a system with the same type of atoms it is simple to generalize it for honeycomb 
lattices with different type of atoms as in BNC hexagonal sheets [24J. However, the 
main differences will be: (i) the number of optical and acoustic branches increases; (ii) 
the different site energies due to different atoms will change the form of the spin wave 
bands. We have also suggested a possible mechanism for ferromagnetism in irradiated 
graphite: the appearance of bound states due to disorder and the enhancement of the 
density of states. If it becomes possible in the future to perform neutron scattering 
experiments on ferromagnetic graphite allotrope the results we present here may have 
direct experimental importance. 
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